To begin with the analysis, we will first load the required packages. The downstream analysis of the scRNA-seq is built upon Seurat and sctransform packages. The scDblFinder package is used for identifying and removing doublets in the dataset. The single-cell annotation is carried out using the marker-based scSorter method.
library("Seurat")
library("ggplot2")
library("sctransform")
library("BiocParallel")
library("scDblFinder")
library("SingleCellExperiment")
library("cowplot")
theme_set(theme_cowplot())
library("RColorBrewer")
library("viridis")
library("scSorter")
library("knitr")
First, we will load the filtered count matrix from cellranger. This matrix contains the scRNA-seq counts from all the HPAP samples.
panc_data <- Read10X(data.dir = "./data/HPAP_samples_aggr_18Apr2022/outs/count/filtered_feature_bc_matrix/")
The dimension of the raw counts matrix after merging all HPAP scRNA-seq samples with cellranger is shown below. There are total of 36601 genes and 269417 cells.
# panc_data <- Read10X(data.dir = "PATH_TO_FEATURE_MATRIX")
dim(panc_data)
## [1] 36601 269417
The Seurat object initialization step above only considered cells that expressed at least 200 genes. The min.cells- parameter will include features detected in at least three cells. The min.features - Include cells where at least three fifty features are detected.
With the above parameters the count matrix will also be subsetted. Optionally, to revive the excluded features, create a new object with a lower threshold.
panc <- CreateSeuratObject(counts = panc_data, min.cells = 3, min.features = 200,
project = "hpap@67")
dim(panc)
## [1] 32066 258379
In the Seurat object, 32066 genes and 258379 cells were present after passing the thresholds.
Here, we will add four new columns to the Seurat metadata. The first column sample_id contains raw sample names from cellranger, hpap_id refers to the donor identifier, the disease_id corresponds to the numeric count in each disease type, and finally the disease_state shows the general state of disease (CTL, AAB, and T1D).
## Add metadata information
sample <- read.csv(file.path("./metadata", "HPAP_samples_aggr_18Apr2022_cellranger.csv"), stringsAsFactors=F)
## Add sample_id
CellsMeta["sample_id"] <- df$sample_id
CellsMetaTrim <- subset(CellsMeta, select = c("sample_id"))
panc <- AddMetaData(panc, CellsMetaTrim)
## Add hpap_id
CellsMeta["hpap_id"] <- df$hpap_id
CellsMetaTrim <- subset(CellsMeta, select = c("hpap_id"))
panc <- AddMetaData(panc, CellsMetaTrim)
## Add disease_id
CellsMeta["disease_id"] <- df$disease_id
CellsMetaTrim <- subset(CellsMeta, select = c("disease_id"))
panc <- AddMetaData(panc, CellsMetaTrim)
## Add disease_state
CellsMeta["disease_state"] <- df$disease_state
CellsMetaTrim <- subset(CellsMeta, select = c("disease_state"))
panc <- AddMetaData(panc, CellsMetaTrim)
head(panc@meta.data) %>%
kable(format = "html", col.names = colnames(head(panc@meta.data))) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "100px")
| orig.ident | nCount_RNA | nFeature_RNA | sample_id | hpap_id | disease_id | disease_state | |
|---|---|---|---|---|---|---|---|
| AAACCTGAGCGTTGCC-1 | hpap@67 | 16875 | 4080 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGAGGACAGAA-1 | hpap@67 | 6497 | 892 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGAGTCGTACT-1 | hpap@67 | 1169 | 253 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGCAACTGCTA-1 | hpap@67 | 13689 | 3553 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGCAATGTAAG-1 | hpap@67 | 16375 | 3583 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGCACCAGGTC-1 | hpap@67 | 16431 | 3438 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
## Calculate MT %
panc[["percent.mt"]] <- PercentageFeatureSet(panc, pattern = "^MT-")
VlnPlot(panc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0.001, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
VlnPlot(panc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
# FeatureScatter used to visualize feature-feature relationships.
plot1 <- FeatureScatter(panc, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "disease_state") & theme(legend.position = 'none')
plot2 <- FeatureScatter(panc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "disease_state")
CombinePlots(plots = list(plot1, plot2))
The scDblFinder function is employed to remove the potential doublets from each HPAP donor separately. The parameter sample_id refers to individual donor. The function uses an object from SingleCellExperiment class, therefore the SingleCellExperiment package is used for conversion.
## Convert object into singlecellexperiment
panc.sce <- as.SingleCellExperiment(panc)
panc.sce <- scDblFinder(panc.sce, samples="sample_id", clusters=FALSE, BPPARAM=MulticoreParam(10))
## Convert sce object back to seurat
panc_seurat <- as.Seurat(panc.sce, counts = "counts", data = "logcounts")
table(panc.sce$scDblFinder.class)
##
## singlet doublet
## 234285 24094
The doublet information is stored in the metadata. The updated metadata can be found below
head(panc@meta.data) %>%
kable(format = "html", col.names = colnames(head(panc@meta.data))) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "100px")
| orig.ident | nCount_RNA | nFeature_RNA | sample_id | hpap_id | disease_id | disease_state | percent.mt | |
|---|---|---|---|---|---|---|---|---|
| AAACCTGAGCGTTGCC-1 | hpap@67 | 16875 | 4080 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 7.330370 |
| AAACCTGAGGACAGAA-1 | hpap@67 | 6497 | 892 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 1.062029 |
| AAACCTGAGTCGTACT-1 | hpap@67 | 1169 | 253 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 4.362703 |
| AAACCTGCAACTGCTA-1 | hpap@67 | 13689 | 3553 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 3.382278 |
| AAACCTGCAATGTAAG-1 | hpap@67 | 16375 | 3583 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 3.877863 |
| AAACCTGCACCAGGTC-1 | hpap@67 | 16431 | 3438 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 1.545858 |
The nFeature_RNA and nCount_RNA represents the number of genes detected in each cell and the total number of molecules (UMIs) detected within a cell respectively. While the low nFeature_RNA in a cell mean that it may be dead/dying or it may represent an empty droplet, the high nCount_RNA and nFeature_RNA denotes that the cell may be a doublet or multiplet. This filtering along with mitochondrial reads is crucial pre-processing step, because, removing such outliers from these groups might also remove some of the doublets or dead/empty droplets.
The cell being either doublet or singlet is stored in column scDlFinder.class. The potential doublets are removed in this step.
# After doublet removal
panc <- subset(panc_seurat, subset = scDblFinder.class == "singlet")
panc_new <- subset(panc, subset = nFeature_RNA > 200 & nFeature_RNA < 9000
& percent.mt < 25
& nCount_RNA < 120000
)
dim(panc_new)
## [1] 32066 222679
After removing doublets, now we can proceed with further pre-processing steps. For ex, Mitochondrial genes are useful indicators of cell state.
Adding the Mitochondrial genes
panc[["percent.mt"]] <- PercentageFeatureSet(panc, pattern = "^MT-")
The below plot shows the number of features, counts, and the percentage of mitochondrial reads after filtering.
VlnPlot(panc_new, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0.001, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
VlnPlot(panc_new, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
The FeatureScatter is typically used to visualize feature-feature relationships.
plot1 <- FeatureScatter(panc_new, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "disease_state") & theme(legend.position = 'none')
plot2 <- FeatureScatter(panc_new, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "disease_state")
CombinePlots(plots = list(plot1, plot2))
The seurat’s SCTranform (SCT) function is used to normalize the counts that measures the differences in sequencing depth per cell for each sample. The SCT method is built on the concept regularized negative binomial model to perform normalization and variance stabilization of single-cell RNA-seq data. It removes the variation due to sequencing depth (nUMIs). The vars.to.regress parameter allows us to regress out the variation from other sources such as percentage of mitochondrial. The output of SCT model is the normalized expression levels for all the transcripts. Lastly, the variable.features.n parameter is used to select the variable features and is set as 3000.
panc_new_sct <- SCTransform(panc_new, method = "glmGamPoi", variable.features.n = 3000, vars.to.regress = "percent.mt", verbose = FALSE)
The seurat object now consists of two assays RNA and SCT, with the default being set as SCT.
DefaultAssay(panc_new_sct)
## [1] "SCT"
These are now standard steps in the Seurat workflow for visualization and clustering
panc_new_sct <- RunPCA(panc_new_sct, verbose = TRUE)
Plot showing heatmaps of PCs- 1 to 12 from the PCA analysis
DimHeatmap(panc_new_sct, dims = 1:12, cells = 500, balanced = TRUE)
Plot showing heatmaps of PCs- 13 to 24 from the PCA analysis
DimHeatmap(panc_new_sct, dims = 13:24, cells = 500, balanced = TRUE)
Plot showing heatmaps of PCs- 25 to 36 from the PCA analysis
DimHeatmap(panc_new_sct, dims = 25:36, cells = 500, balanced = TRUE)
Plot showing the variance explained by PCs from the PCA analysis
ElbowPlot(panc_new_sct, ndims = 40)
ElbowPlot(panc_new_sct)
The Uniform Manifold Approximation and Projection (UMAP) dimensional reduction technique is run on the first 20 principal components.
panc_new_sct <- RunUMAP(panc_new_sct, dims = 1:20, verbose = FALSE)
Here we compute the k nearest neighbors in the dataset. The nearest neighbors are calculated on the reduced dimensions from the PCA. We use first 20 principal components.
panc_new_sct <- FindNeighbors(panc_new_sct, dims = 1:20, verbose = FALSE)
Once we have calculated the k-nearest neighbors in previous step, we will find the clusters of cells by the shared nearest neihbor (SNN) clustering algorithm. Here, we have set the resolution as 0.8. The higher resolution usually gives more number of clusters.
panc_new_sct <- FindClusters(panc_new_sct, resolution = 0.8, method = "igraph", verbose = FALSE)
UMAP plot showing the clusters based on SCT data and snn optimization at resolution 0.8
DimPlot(panc_new_sct, group.by='SCT_snn_res.0.8', reduction='umap') +
ggtitle('SCT_snn_res.0.8')
UMAP plot showing the distribution of cells across different disease types.
DimPlot(panc_new_sct, group.by='disease_state', reduction='umap') +
ggtitle('Disease type')
UMAP plot showing the cell distribution with respect to individual donors
DimPlot(panc_new_sct, group.by='hpap_id', reduction='umap') +
ggtitle('Sample type')
Loading marker genes data
anno <- read.csv("./annot/anno.csv")
anno[,1] <- NULL
We used marker-based scSorter cell annotation method for assigning cells to known cell type based on the marker genes. The markers for each cell type were added based on comprehensive literature search.
anno %>%
kable(format = "html", col.names = colnames(anno)) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "100px")
| Type | Marker |
|---|---|
| Acinar | PRSS1 |
| Acinar | REG1A |
| Acinar | CPA1 |
| Acinar | CPA2 |
| Alpha | GCG |
| Alpha | GC |
| Alpha | TTR |
| Beta | INS |
| Beta | IGF2 |
| Beta | IAPP |
| Ductal | KRT19 |
| Ductal | CFTR |
| Ductal | SFRP5 |
| Ductal | MMP7 |
| Stellates_Mesenchymal | COL1A1 |
| Stellates_Mesenchymal | PDGFRB |
| Stellates_Mesenchymal | RGS10 |
| Stellates_Mesenchymal | THY1 |
| Endothelial | VWF |
| Endothelial | CD93 |
| Immune | NCF2 |
| Immune | PTPRC |
| PP_Gamma | PPY |
| PP_Gamma | CARTPT |
| PP_Gamma | PCDH10 |
| PP_Gamma | PLAC8 |
| Delta | SST |
| Delta | LEPR |
| Delta | PRG4 |
| Epsilon | GHRL |
The scSorter method is based on a semi-supervised learning algorithm. It expects the pre-processed data as input. Therefore, we used the previously normalized, scaled and transformed expression matrix generated by SCT method as input. The top 3000 highly variable genes selected by SCT method was also given as input for scSorter method.
expr <- panc_new_sct@assays$SCT@data
topgenes <- head(panc_new_sct@assays$SCT@var.features, 3000)
topgene_filter = rowSums(expr[topgenes, ]!=0) > ncol(expr)*.1
topgenes = topgenes[topgene_filter]
## At last, we subset the preprocessed expression data and run scSorter.
picked_genes = unique(c(anno$Marker, topgenes))
expr = expr[rownames(expr) %in% picked_genes, ]
Based on the expression matrix and the marker genes the scSorter method is fit and the cell type assignment results predicted are stored in the Pred_Type vector.
rts <- scSorter(expr, anno)
panc_new_sct$cell_type <- rts$Pred_Type
table(rts$Pred_Type)
##
## Acinar Alpha Beta
## 59496 65075 45009
## Delta Ductal Endothelial
## 5664 20300 5864
## Epsilon Immune PP_Gamma
## 196 1220 2821
## Stellates_Mesenchymal Unknown
## 9583 7451
The celltype information is added to the metadata in the original Seurat object
cell_type.info <- data.frame(cell_type = panc_new_sct$cell_type, row.names= colnames(panc_new_sct_umap))
panc_new_sct_umap <- AddMetaData(object = panc_new_sct_umap, metadata = cell_type.info)
## Convert cell type to factor
panc_new_sct_umap$cell_type <- as.factor(panc_new_sct_umap$cell_type)
Idents(panc_new_sct_umap) <- "cell_type"
panc_new_sct_umap$ident <- NULL
FeaturePlot(panc_new_sct_umap, features = c("PRSS1", "REG1A", "CPA1", "CPA2", "GCG",
"GC", "TTR", "INS", "IAPP", "GHRL",
"COL1A1", "PDGFRB", "RGS10", "THY1", "VWF",
"CD93", "PPY", "SST", "NCF2", "PTPRC"), pt.size = 0.2, ncol = 5) + theme(legend.position = 'none',
axis.title.x = element_blank(), axis.title.y = element_blank())
The UMAP showing the cell type annotations from the scSorter method
DimPlot(panc_new_sct_umap, group.by = "cell_type", reduction = "umap", label = FALSE, cols= c("#e30800", "#f56505", "#dec400", "#006630", "#0223c7","#5b02c7", "#00b0e6", "#c40080", "#02f00a", "#7d3301", "#000000"))
UMAP showing only Acinar cells
DimPlot(panc_new_sct_umap, group.by = "cell_type", reduction = "umap", label = FALSE, cols= c("#fac720", "#000000", "#000000", "#000000", "#000000","#000000", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000"))
UMAP showing only Alpha cells
UMAP showing only Beta cells
UMAP showing only Ductal cells
UMAP showing only Immune cells